$\eta = 0.008$
# for post jump HJB
import os
import numpy as np
import pandas as pd
import pickle
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import plotly.graph_objects as go
mpl.rcParams["savefig.bbox"] = "tight"
with open("data/res-28-14-24", "rb") as file:
data = pickle.load(file)
Kd = data["Kd"]
Kd_max = data["Kd_max"]
Kg = data["Kg"]
Kg_max = data["Kg_max"]
Kd_mat = data["Kd_mat"]
Kd_max = data["Kd_max"]
Kg_mat = data["Kg_mat"]
Kg_max = data["Kg_max"]
Y = data["Y"]
Y_mat = data["Y_mat"]
i_d = data["i_d"]
i_g = data["i_g"]
v = data["v0"]
print(v.shape)
(36, 99, 41)
fig = go.Figure(data=[go.Surface(z=v[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0] ),
go.Surface(z=v[:,:,25],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showscale=False),
go.Surface(z=v[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showscale=False),
])
fig.update_layout(title='V(Kd, Kg), for given Y',
width=1000, height=800,
margin=dict(l=65, r=50, b=10, t=90),
scene = dict(
xaxis_title='Total capital',
yaxis_title='Green capital / total capital ratio',
zaxis_title='Value function'),
)
fig.show()
fig = go.Figure(data=[go.Surface(z=i_d[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0],
showlegend=True,
name="Y=0",
showscale=False ),
go.Surface(z=i_d[:,:,20],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showlegend=True,
name="Y=2",showscale=False),
go.Surface(z=i_d[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], name="Y=4",
showlegend=True,
showscale=False),
])
fig.update_layout(title='i_d(K, r), for given Y',
width=1000, height=800,
margin=dict(l=65, r=50, b=10, t=90),
legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01),
scene = dict(
xaxis_title='Total capital',
yaxis_title='Green capital / total capital ratio',
zaxis_title='Dirty investment'),
)
fig.show()
fig = go.Figure(data=[go.Surface(z=i_g[:,:,0], x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0],
showlegend=True,
name="Y=0",
showscale=False ),
go.Surface(z=i_g[:,:,20],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], showlegend=True,
name="Y=2",showscale=False),
go.Surface(z=i_g[:,:,-1],x=np.exp(Kd_mat[:,:,0]), y=Kg_mat[:,:,0], name="Y=4",
showlegend=True,
showscale=False),
])
fig.update_layout(title='i_g(K, r), for given Y',
width=1000, height=800,
margin=dict(l=65, r=50, b=10, t=90),
legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01),
scene = dict(
xaxis_title='Total capital',
yaxis_title='Green capital / total capital ratio',
zaxis_title='Green investment'),
)
fig.show()
plt.plot(np.exp(Kd )* Kg[2], v[:, 2,25], label="Green capital / Total capital = ".format(Kg[2]))
plt.plot(np.exp(Kd) * Kg[25], v[:, 25,25], label="Green capital / Total capital = ".format(Kg[25]))
plt.plot(np.exp(Kd) * Kg[-1], v[:, -1,25], label="Green capital / Total capital = ".format(Kg[-1]))
plt.title("Value function")
plt.xlabel("Green capital")
Text(0.5, 0, 'Green capital')
plt.plot(np.exp(Kd), i_d[:, 0,20], label="Clean capital / Total capital = {}".format(Kg[0]) )
plt.plot(np.exp(Kd), i_d[:, 25,20], label="Clean capital / Total capital = {}".format(Kg[25]))
plt.plot(np.exp(Kd), i_d[:, -1,20], label="Clean capital / Total capital = {}".format(Kg[-1]))
plt.ylim(0)
plt.xlabel("Total capital")
plt.title("Dirty investment, Y = 2")
plt.legend()
<matplotlib.legend.Legend at 0x7f9dabe98be0>
plt.plot(np.exp(Kd), i_d[:, 0,40], label="Clean capital / Total capital = {}".format(Kg[0]))
plt.plot(np.exp(Kd), i_d[:, 25,40], label="Clean capital / Total capital = {}".format(Kg[25]))
plt.plot(np.exp(Kd), i_d[:, -1,40], label="Clean capital / Total capital = {}".format(Kg[-1]))
plt.xlabel("Total capital")
plt.title("Dirty investment, Y = 4")
plt.ylim(0)
plt.legend()
<matplotlib.legend.Legend at 0x7f9dabcfd310>
plt.plot(Y, v[25, 0, :], label="Green capital / Total capital = ".format(Kd[0]))
plt.plot(Y, v[25, 25,:], label="Green capital / Total capital = ".format(Kd[25]))
plt.plot(Y, v[25, -1,:], label="Green capital / Total capital = ".format(Kd[-1]))
# plt.ylim(0)
plt.title("Value function")
plt.xlabel("Temperature anomaly")
Text(0.5, 0, 'Temperature anomaly')
plt.plot(Y, i_d[0, 25, :])
plt.plot(Y, i_d[25, 25,:])
plt.plot(Y, i_d[-1, 25,:])
plt.title("Dirty investment")
plt.xlabel("Temperature anomaly")
plt.ylim(0)
(0.0, 0.08542750080854612)
# fig 1
plt.figure()
plt.title("value function - total capital")
plt.plot(np.exp(Kd), v[:, 10,10], label="Green capital / Total capital = {:.2f}, $Y = {:.2f}$".format(Kg[10], Y[10]))
plt.xlabel("Total capital")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Kd.pdf")
plt.show()
# fig 2
plt.figure()
plt.title("value function - green capital")
plt.plot(np.exp(Kd[30]) * Kg, v[30,:,20], label="Total capital = {:d}, $Y = {:.2f}$".format(int(np.exp(Kd[30])), Y[20]))
plt.xlabel("Green capital")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Kg.pdf")
plt.show()
# fig 3
plt.figure()
plt.title("value function - temperature")
plt.plot(Y, v[20,20,:], label="Total capital = {:d},\n green capital / total capital = {:.2f}".format(int(np.exp(Kd[20])), Kg[20]))
# plt.plot(Y, v[5,5,:], label="$K_d = {:d}, K_g = {:d}$".format(int(Kd[5]), int(Kg[5])))
plt.xlabel("$Y$")
plt.ylabel("v")
plt.legend()
# plt.savefig(figuredir + "v-Y.pdf")
plt.show()
# fig 1
plt.figure()
plt.title("dirty investment - total capital")
plt.plot(np.exp(Kd), i_d[:,50,20],
label="green capital / total capital = {:.2f}, $Y = {:.2f}$".format(Kg[50], Y[20]))
plt.xlabel("Total capital")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Kd.pdf")
plt.show()
# fig 2
plt.figure()
plt.title("dirty investment - green capital / total capital")
plt.plot(Kg, i_d[20,:,20], label="$K_d = {:d}, Y = {:.2f}$".format(int(np.exp(Kd[20])), Y[20]))
plt.xlabel("green capital / total capital")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Kg.pdf")
plt.show()
# fig 3
plt.figure()
plt.title("dirty investment - temperature")
plt.plot(Y, i_d[20,20,:], label="$K_d = {:d}, K_g = {:d}$".format(int(np.exp(Kd[20])), int(Kg[20])))
plt.xlabel("$Y$")
plt.ylabel("$i_d$")
plt.legend()
# plt.savefig(figuredir + "id-Y.pdf")
plt.show()
# print(i_g)
# fig 1
plt.figure()
plt.title("green investment - total capital")
plt.plot(np.exp(Kd), i_g[:,50,20], label="green capital / total = {:.2f}, $Y = {:.2f}$".format(Kg[50], Y[20]))
plt.xlabel("total capital")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Kd.pdf")
plt.show()
# fig 2
plt.figure()
plt.title("green investment - total capital")
plt.plot(Kg, i_g[20,:,20], label="total capital = {:d}, $Y = {:.2f}$".format(int(np.exp(Kd[20])), Y[20]))
plt.xlabel("green capital / total capital")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Kg.pdf")
plt.show()
# fig 3
plt.figure()
plt.title("green investment - temperature")
plt.plot(Y, i_g[20,50,:],
label="total capital = {:d}, \n green capital / total capital = {:.2f}".format(int(Kd[20]), Kg[50]))
plt.xlabel("$Y$")
plt.ylabel("$i_g$")
plt.legend()
# plt.savefig(figuredir + "ig-Y.pdf")
plt.show()